library(BiodiversityR)
library(raster)
library(rgdal)
library(dismo)
library(terra)
library(ggmap)
library(ggspatial)
library(ggforce)

files = list.files("data/grds", pattern = "grd", 
                   full.names = TRUE)
predictors.raster = stack(files)

predictors.raster = raster::subset(predictors.raster, 
                                   subset=c("bio5", "bio6", "bio16", "bio17"))
predictors.raster

VIF.recheck = ensemble.VIF(predictors.raster)


files = list.files("data/tifs", pattern = "tif", 
                   full.names = TRUE)
predictors = rast(files) %>% 
  terra::subset(subset = c("bio5","bio6","bio16","bio17"))

predictors  # 同predictors.raster
plot(predictors)

pres = read.csv("data/bradypus.csv")[-1]
pres

pres.spat = ensemble.spatialThin(x = pres, thin.km = 10)
pres.env = ensemble.environmentalThin(x = pres, 
                                      predictors.stack = predictors.raster,
                                      thin.n = nrow(pres.spat))   # 108个

# environmentally-thinned数据集，最小环境距离为0
# 物种点进一步减少到90
pres.env = ensemble.environmentalThin(x = pres,
                                      predictors.stack = predictors.raster,
                                      thin.n = 90)

bbox_study2 = c(left = min(pres$lon - 5),
                bottom = min(pres$lat - 5),
                right = max(pres$lon + 8),
                top = max(pres$lat + 5))
stamen_study2 = get_stamenmap(bbox_study2, zoom=5, force=TRUE,
                              maptype="terrain-background")
stamen_map2 = ggmap(stamen_study2)
map_study2 = stamen_map2 +
  geom_mark_hull(data = pres,
                 mapping = aes(x=lon, y=lat),
                 colour = "red") +
  geom_point(data = pres.env,
             aes(x=lon, y=lat), colour="black", 
             size=2.0, shape=21, fill="green") +
  geom_point(data = dplyr::setdiff(pres, pres.env),
             aes(x=lon, y=lat), colour="red", 
             alpha=0.5, size=3.0, shape=4, stroke=2) +
  theme(axis.title = element_blank())
map_study2

circles.model = circles(p = pres.env, d = 500000,
                        lonlat = TRUE)
circles.predicted = predict(circles.model,
                            predictors.raster[[1]])
circles.data = as.data.frame(circles.predicted,
                             xy = TRUE, na.rm = TRUE)
map_study2 = stamen_map2 +
  geom_contour_filled(data = circles.data,
                      aes(x=x, y=y, z=layer),
                      colour="black", alpha=0.4,
                      show.legend=FALSE) +
  geom_mark_hull(data=pres,
                 mapping = aes(x=lon, y=lat),
                 colour="red") +
  geom_point(data=pres.env,
             aes(x=lon, y=lat), colour="black", 
             size = 2.0, shape = 21, fill = "red") +
  theme(axis.title = element_blank())
map_study2

background1 <- dismo::randomPoints(mask=circles.predicted,
                                   n=2000,
                                   p=pres.env,
                                   excludep=TRUE)
background1.data = terra::extract(predictors, 
                                  y = background1)
background2 = background1[complete.cases(background1.data),]
background3 = background2[sample(nrow(background2), 
                                 size = 10*nrow(pres.env)),]
background3 = data.frame(background3)
names(background3) = c("lon", "lat")

locations.folded = ensemble.spatialBlock(
  x = predictors.raster, p = pres.env,
  a = background3, theRange = 500000,
  return.object = TRUE)

pres.folded = data.frame(pres.env,
                         fold=as.character(locations.folded$k$groupp))
background.folded = data.frame(background3,
                               fold=as.character(locations.folded$k$groupa))

map_study2 = stamen_map2 +
  geom_mark_hull(data=pres,
                 mapping = aes(x=lon, y=lat),
                 colour="red") +
  geom_point(data=pres.folded,
             aes(x=lon, y=lat, fill=fold),
             colour="black", size=2, shape=21) +
  theme(axis.title = element_blank())
map_study2

map_study2 = stamen_map2 +
  geom_mark_hull(data=pres,
                 mapping = aes(x=lon, y=lat),
                 colour="red") +
  geom_point(data=background.folded,
             aes(x=lon, y=lat, fill=fold),
             colour="black", size=1.2, shape=21) +
  theme(axis.title = element_blank())
map_study2

ensemble.calibrate.step1 = ensemble.calibrate.weights(
  x = predictors, p = pres.env, a = background3,
  k = locations.folded$k,
  SINK=FALSE, species.name = "Bradypus",
  MAXENT=0, MAXNET=1, MAXLIKE=1, GBM=0, GBMSTEP=1, 
  RF=1, CF=0, GLM=1, GLMSTEP=1, GAM=1, GAMSTEP=1, 
  MGCV=1, MGCVFIX=0, EARTH=1, RPART=1, NNET=1, 
  FDA=1, SVM=1, SVME=1, GLMNET=1, BIOCLIM.O=0, 
  BIOCLIM=1, DOMAIN=1, MAHAL=1, MAHAL01=0,
  ENSEMBLE.tune=TRUE, PROBIT=TRUE,
  ENSEMBLE.best=0, ENSEMBLE.exponent=c(1, 1.5, 2),
  ENSEMBLE.min=0.55,
  Yweights="BIOMOD",
  PLOTS=FALSE, formulae.defaults=TRUE)

model.weights = ensemble.calibrate.step1$output.weights
model.weights
x.batch = ensemble.calibrate.step1$x
p.batch = ensemble.calibrate.step1$p
a.batch = ensemble.calibrate.step1$a

ensemble.calibrate.step2 = ensemble.calibrate.models(
  x=x.batch, p=p.batch, a=a.batch,
  SINK=FALSE, species.name="Bradypus",
  models.keep=TRUE,
  input.weights=model.weights,
  ENSEMBLE.tune=FALSE, PROBIT=TRUE,
  Yweights="BIOMOD",
  PLOTS=FALSE, formulae.defaults=TRUE)

ensemble.terra.results = ensemble.terra(
  xn=predictors.terra,                                             models.list=ensemble.calibrate.step2$models,                     input.weights=model.weights,
  SINK=FALSE, evaluate=TRUE,                                       RASTER.species.name="Bradypus_terra", 
  RASTER.stack.name="base")

suit.file = "ensembles/suitability/Bradypus_terra_base.tif"
suit.raster = terra::rast(suit.file) / 1000
plot(suit.raster)

suit.raster.masked = terra::mask(suit.raster,
                                 mask=pres.raster,
                                 maskvalue=0,
                                 updatevalue=NA) 
# 计算分位数
breaks1 = quantile(as.data.frame(suit.raster.masked,
                                 xy=FALSE, na.rm=TRUE)[,1])
p_values = terra::extract(suit.raster, y=pres.env)[,2]
a_values = terra::extract(suit.raster, y=background3)[,2]

breaks2 = dismo::threshold(
  dismo::evaluate(p=p_values, a=a_values),
  stat = "sensitivity", sensitivity=0.90)
breaks1 = c(0.0, breaks2, breaks1)
suit.classified = terra::classify(suit.raster, 
                                  rcl=as.matrix(breaks1),
                                  include.lowest=TRUE)
plot.colours = RColorBrewer::brewer.pal(
  length(breaks1), "Greens")
plot.colours[1:2] = c("grey60","#D8B365")
plot(suit.classified, col=plot.colours)
